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Abstract. 

We  study  how  the  Newton-GMRES  iteration  can  enable  dynamic  simulators  (time-steppers)  to 
perform  fixed-point  and  path-following  computations.  For  a  class  of  dissipative  problems,  whose 
dynamics  are  characterized  by  a  slow  manifold,  the  Jacobian  matrices  in  such  computations  are 
compact  perturbations  of  the  identity.  We  examine  the  number  of  GMRES  iterations  required  for 
each  nonlinear  iteration  as  a  function  of  the  dimension  of  the  slow  subspace  and  the  time-stepper 
reporting  horizon.  In  a  path-following  computation,  only  a  small  number  (one  or  two)  of  additional 
GMRES  iterations  is  required. 
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1.  Introduction.  In  studying  the  dynamic  behavior  of  evolution  equations, 
(1.1)  du/dt  =  f(u;  X), 

a  computational  modeler  typically  chooses  between  two  paths:  the  first  is  developing 
a  dynamic  simulator  for  the  process;  the  second  is  developing  algorithms  to  locate 
particular  features  of  the  long-term  dynamics  of  the  process,  such  as  steady  states 
or  limit  cycles.  The  first  path  typically  gives  rise  to  initial  value  problems,  and  the 
corresponding  codes  are  dynamic  simulators  -  we  will  call  them  “time-steppers”  since, 
given  the  state  of  the  system  at  a  moment  in  time  they  produce  (an  approximation 
of)  the  state  of  the  system  at  a  later  moment.  The  second  path  typically  gives  rise  to 
fixed  point  algorithms  for  solving  coupled  nonlinear  algebraic  equations;  these  may 
be  the  steady  state  equations  themselves,  or  an  augmented  set  arising  in  a  continua¬ 
tion/bifurcation  context.  The  Recursive  Projection  Method  of  Shroff  and  Keller  [23] 
(see  also  the  Adaptive  Condensation  of  Jarausch  and  Mackens  [12,13]  for  symmetric 
problems)  is  an  example  of  an  algorithm  that,  in  some  sense,  connects  the  two  paths. 
The  main  idea  is  to  construct  a  computational  superstructure  that  designs  and  com¬ 
bines  several  calls  to  an  existing  (“legacy”)  time-stepper,  effectively  turning  it  into  a 
fixed  point  solver  for 


u  —  A)  =  0, 

where  <I>t  is  the  result  of  integration  of  (1.1)  with  initial  condition  u  for  time  T  (see 
the  review  of  Tuckerman  and  Barkley  [29]).  Here  the  action  of  the  linearization  of  the 
time-stepper  is  estimated  in  a  matrix-free  fashion  by  the  integration  of  appropriately 
chosen  nearby  initial  conditions. 

Over  the  last  few  years,  this  matrix-free  computational  enabling  technology  has 
found  new  applications  in  the  field  of  multiscale  computations.  In  current  modelling 
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practice,  dynamic  models  are  often  constructed  at  a  microscopic/stochastic  level  of  de¬ 
scription  (e.g.  molecular  dynamics,  kinetic  Monte  Carlo  or  Lattice  Boltzmann  codes); 
the  closure  required  to  obtain  explicit,  macroscopic  “effective”  equations  is  not  avail¬ 
able  in  closed  form.  Algorithms  like  RPM  can  then  be  “wrapped  around”  appropri¬ 
ately  initialized  ensembles  of  short  bursts  of  microscopic  simulations,  and  solve  the 
effective  equation  without  ever  obtaining  it  in  closed  form.  Several  examples  of  such 
“coarse  equation-free”  computation  have  now  been  explored,  and  the  mathematical 
underpinnings  of  the  approach  are  being  extensively  studied  [9,19,25,27]. 

In  this  paper  we  study  the  Newton-GMRES  iteration  as  a  computational  “wrap¬ 
per”  around  a  legacy  time-stepper.  This  wrapper  enables  the  computation  and  contin¬ 
uation  of  fixed  points  of  the  time-T  map  of  the  time-stepper  (i.e.  steady  states  of  the 
corresponding  dynamical  equations).  It  is  useful,  for  purposes  of  discussion,  to  con¬ 
sider  that  the  time-stepper  is  available  as  an  input-output  black  box  (an  executable) 
which  cannot  be  modified. 

The  paper  is  organized  as  follows:  We  first  review  certain  properties  of  GMRES 
in  the  context  of  problems  whose  linearization  is  a  compact  perturbation  of  the  iden¬ 
tity.  We  argue  that  time-steppers  for  a  class  of  dissipative  problems,  whose  long-term 
dynamics  lie  on  a  low-dimensional,  slow  manifold,  may  fit  this  description;  we  then 
proceed  to  examine  Newton-GMRES  convergence  and  number  of  iterations  for  fixed 
point  computation  of  such  time-steppers  in  a  continuation  context.  Our  numerical  ex¬ 
amples  are  the  Chandrasekhar  H-equation  [5]  and  the  time-stepper  of  a  discretization 
of  a  reaction-diffusion  problem  known  to  possess  a  low-dimensional  inertial  mani¬ 
fold  [8,26].  We  illustrate  the  effect,  on  the  GMRES  iterations,  of  the  time-stepper 
reporting  horizon  both  for  fixed  point  and  for  pseudo-arclength  continuation  compu¬ 
tations.  We  conclude  with  a  brief  discussion  and  thoughts  about  extensions  of  this 
approach. 

2.  Convergence  Analysis. 

2.1.  GMRES  Preliminaries.  We  will  use  Newton-GMRES  to  solve  the  non¬ 
linear  fixed  point  equation 

(2.1)  u  —  $t(w;  A)  =  F(u)  =  0. 

In  what  follows  we  will  work  in  RN  and  use  the  usual  Euclidean  norm;  the  idea  is  that 
equation  (2.1)  arises  from  integrating  a  set  of  ODEs,  possibly  from  the  discretization 
of  a  partial  differential  equation  on  a  given  mesh.  We  will  also  assume  (and  the 
consequences  of  this  will  become  apparent  below)  that  the  long-term  dynamics  of  the 
discretized  PDE  occur  on  an  attracting  “slow  manifold”  of  low  dimension  p  «  N-, 
p  will  be  assumed  independent  of  mesh  refinement  (i.e.  it  will  remain  constant  as  N 
increases.) 

GMRES  is  an  iterative  method  for  solving  linear  systems  Ax  =  b  in  RN .  The 
kth  GMRES  iteration  minimizes  the  residual  r  over  Xq  +  /Cfc,  where  Xq  is  the  initial 
iterate  and  ICj-  is  the  the  fcth  Krylov  subspace 

/Cfc  =  span{r0,  Ar0, ...,  Ak~lr0}. 

A  consequence  [15, 22]  of  the  minimization  is  that 

(2.2)  ||rfc||  =  min  ||p(A)r0|| 

p&vk 

where  Vk  is  the  space  of  kth  degree  residual  polynomials,  i.e.  polynomials  of  degree 
k  such  that  p(0)  =  1. 
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We  will  apply  (2.2)  to  a  special  class  of  problems  where 
(2.3)  A  =  I  -  K  +  E. 

In  (2.3) 

•  I  —  K  is  nonsingular, 

•  K  =  PdKPd,  where  Pd  is  an  orthogonal  projection  onto  a  space  D  of 
dimension  p  «  N,  and 

•  E  is  a  matrix  with  small  norm. 

We  will  analyze  the  performance  of  GMRES  in  a  way  different  from  the  eigenvalue- 
based  approach  used  for  cliagonalizable  matrices  [10, 15, 18,28].  For  the  class  of  prob¬ 
lems  of  interest  here,  we  can  prove  a  convergence  result  directly  using  methods  similar 
to  those  in  [3,4,7,17].  The  result  in  this  paper  has  stronger  and  sharper  convergence 
rates.  This  result  carries  through  in  the  infinite-dimensional  case  also,  using  the  L2 
norm  for  the  corresponding  function  spaces. 

Theorem  2.1.  Let  A  be  given  by  (2.3).  Then  there  exists  C  such  that  for  all 
m  >  1, 

(2-4)  \\rm{p+1)\\<Cm\\E\r. 

Proof.  Let  pc  be  the  characteristic  polynomial  of  I  —  K.  Since  D  has  dimension 
p,  pc  has  degree  p  +  1.  Clearly 

Pc(A)  =  pc {I  -  K)  +  A  =  A. 

by  the  Cayley-Hamilton  theorem.  Moreover,  there  is  C  >  0  such  that 

l|A||<C'||E||. 

In  fact,  if 

p+  i 

p{z)  =  \  +  Y^lkZk 
k= 1 

then  we  may  use 

p+i 

C=J2  -  PdKPdW*-1  +  0(\\E\\2). 

k= 1 

Hence 

\\p(A)\\<C\\E\\. 

This  is  (2.4)  for  m  =  1. 

Define 

P(z )  =Pc{z)/pc(  0). 

Clearly  p  e  Pp+i ■  Hence,  by  (2.2) 

(2-5)  IIMp+dH  <  \\Pm(A)r0\\  <  Cm\\E\n\r0l 
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as  asserted.  □ 

The  estimate  in  (2.5)  does  not  depend  on  the  eigenvalues  of  A,  nor  is  there  any 
requirement  that  A  be  normal  or  even  diagonalizable. 

As  is  standard,  a  GMRES  iteration  is  terminated  when 

lkfc||  <  r?||r0||, 

where  77  is  an  user-defined  parameter.  In  the  context  of  this  paper  ||IS||  is  well  below 
the  termination  tolerance  r,  so  we  conclude  that  the  iteration  will  terminate  in  at 
most  p  +  1  iterations.  In  the  general  case,  of  course,  more  cycles  could  be  required. 

2.2.  Time-Steppers  and  Steady  State  Solutions.  Let  <f>(u)  denote  the  out¬ 
put  of  the  time-stepper  with  time  step  T  and  initial  data  u:  we  have  dropped  the 
subscript  T  of  $  for  convenience.  We  seek  to  solve 

(2.6)  F(u)  =  u-  $(u)  =  0, 

to  find  a  steady  state  solution.  Consider  the  structure  of  the  eigenvalues  pt  =  exp (cuT) 
of  <f>u  given  the  structure  of  the  eigenvalues  cq  of  the  linearization  of  the  original 
problem.  We  assume  that  there  exists  a  significant  gap  between  a  few  ct,;  close  to  zero 
( p  of  them,  to  be  exact)  and  the  remaining  large  negative  ay .  The  p  small  norm  er,  do 
not  have  to  be  stable  -  they  can  also  be  slightly  unstable,  so  that  the  corresponding 
Pi  are  close  to  the  unit  circle. 

We  then  assume  that  T  is  large  enough  so  that  there  is  a  space  D  of  low  dimension 
p,  so  that 

$(u)  =  Pd^(Pdu)  +  E(u) 

and  E  and  its  Jacobian  Eu  are  small.  For  a  problem  where  some  of  the  p  small  a, 
are  unstable,  this  large  time  T  should  not  be  so  large  that  the  unstable  modes  will 
cause  the  trajectory  to  move  significantly  away  from  the  fixed  point.  In  such  cases 
the  proper  choice  of  T  is  a  somewhat  delicate  matter  [19,23]. 

The  equation  for  the  Newton  step  from  a  current  iterate  uc  is 

(2.7)  Fu(uc)s  =  s  -  $u(uc)s  =  -F(uc), 

here  Fu  and  <I>U  are  the  Jacobians  of  F  and  $. 

A  Newton-GMRES  [1,2,15]  method  solves  (2.7)  with  GMRES,  terminating  the 
linear  (or  inner)  iteration  when 

ll-F«K)s||  <  Vc\\F(uc)\\ 

where  rjc  may  be  changed  as  the  outer  (or  nonlinear)  iteration  progresses  [6,15,16]. 
For  the  application  here,  we  assume  that  ||E,(u)||  is  much  less  than  any  choice  of  rj 
we  make  during  the  iteration. 

The  linear  system  (2.7)  fits  exactly  in  to  the  paradigm  of  §  2.1  with  K  = 
Pd$u{Pdu)Pd  and  E  =  Eu(u).  Hence  we  conclude  that  a  GMRES  iteration  for 
(2.7)  will  take  at  most  p  +  1  iterations  to  drive  the  residual  to  0(||£7|| ). 

2.3.  Time-Steppers  and  Continuation.  In  a  continuation  context  <I>  depends 
on  a  parameter  A.  As  is  standard  [14]  we  add  an  additional  arclength  parameter  s  to 
obtain  the  augmented  system 


(2.8) 


G(u,  A,  s ) 


F{u(s),X{s)) 

uT (u  -  u0)  +  A(A  -  A0)  -  (s  -  So) 
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In  (2.8),  u  and  A  are  approximations  to  du/ds  and  dX/ds ,  which  can  be  obtained  in 
several  ways  [14] .  In  the  calculations  reported  in  §  3  we  used  the  estimate  of  the  slope 
given  by  the  last  two  points,  (uq,  Ao)  and  (w_i,  A_i)  computed  on  the  branch. 

G  is  defined  on  RN+1  and 


Gu,x(u,X)  = 


Fu  F\ 
iiT  A 


We  seek  to  show  that  G'  also  fits  into  our  paradigm,  with  p  replaced  by  at  most 
p  +  2.  Hence,  at  most  two  additional  linear  iterations  will  be  needed  for  each  Newton 
iteration  of  the  augmented  system. 

We  use  (2.6)  to  obtain 


Gu,  x  = 


I  —  Pd$u(Pdu-,  X)Pd  —Pd&\ 
uT  A 


—Eu  — Ex 
0  0 


Let 


where, 


and 


A  = 


I  0 
0  1 


-K  +  £ 


1C  = 


Pd^uPd  Pd®\ 
— uT  1  —  A 


£  = 


—Eu  — Ex 
0  0 


D 

0 


This  fits  the  paradigm  of  §  2.1.  To  see  that,  let 
V  =  span 

Clearly  the  range  of  1C 

R(IC)  =  span 


D 

0 


C  R 


cv. 


JV+1 


To  apply  the  results  from  the  previous  section,  however,  we  need  1C  =  PdICPd ,  where 
P-p  is  the  orthogonal  projector  onto  D.  This  is  why  the  extra  dimension  (ft,  0)T  is 
required. 

To  see  this,  let  y  =  ( u ,  p)  C  RN+1  be  orthogonal  to  V.  This  means  that  y  =  (w,  0), 
where  io  is  orthogonal  to  both  D  and  it.  Clearly 


ICy  = 


Pd&u,Pdw 
•  T 

—U  OJ 


so  /C(/  —  P-d)  =  0.  Since  R{K,)  C  V,  we  have  PpIC  =  1C.  Summarizing 


AC  =  PpICPp 
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and 


dim(X>)  <  p  +  2. 

The  dimension  of  T>  can  be  taken  p  +  1  if  u  is  nearly  in  the  range  of  Pd,  i-  e. 
(2-9)  ||(7-  Pd)u||  =  0{\\EU\\  +  ||EA||). 

In  this  case,  we  can  let 


V  =  span 


D 

0 


C  R 


N+l 


and 


/  Pd^uPd  Pd&\  \ 

\  —{Pdu)t  1  -  A  )  ’ 


(  -Eu  -Ex  \ 

\{{I-PD)u)T  0  )■ 


If  Fu  is  well  conditioned  and  |A|  and  ||FA||  are  not  too  large,  then  (2.9)  holds.  To 
see  this  we  differentiate  F  =  0  and  obtain 


Fuu  +  F\X  —  0, 


which  implies,  if  Fu  is  nonsingular,  that 

(2.10)  u  =  —F~1F\\. 

Our  assumptions  are  that  ||£?u||  and  ||S^||  are  much  smaller  than  ||FU||  and  ||Ta||. 
Hence,  if  Fu  is  well  conditioned,  the  Banach  lemma  implies  that 

F-1  =  (J  -  Pd^uPd)-1  +  0(||£;u||). 


Moreover, 


Ta  =  -pd$a  +  o(||7;a||). 

We  incorporate  this  into  (2.10)  to  obtain 

(2.11)  u  =  (I  -  PdQM-'PdQx A  +  O  (|A|(||77U||||7\||  +  ||£;A||))  , 

which  implies  (2.9)  if  |A|  and  ||FA||  are  0(1). 

Summarizing,  if  ||T1U||  and  ||7iA||  are  much  smaller  that  ||F„||  and  ||FAj|  respec¬ 
tively,  |A|  and  ||FA||  are  0(1),  and  if  Fu  is  well  conditioned,  then  (2.9)  holds.  Near 
folds  and  bifurcations,  Fu  is  singular,  and  in  those  cases  we  may  need  to  require  that 
the  dimension  of  V  be  p  +  2. 

3.  Numerical  Results.  All  the  computations  were  done  using  MATLAB  ver¬ 
sion  6.0  Release  12  on  a  PC  with  a  Pentium4  2.53GHz  CPU. 
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3.1.  The  H-equation.  As  a  first  example,  we  solved  a  problem  which  does  not 
arise  in  a  time-stepper  context,  but  for  which  the  Jacobian  fits  our  paradigm  [21]. 
The  solution  and  continuation  problem  for  a  100-node  midpoint  rule  discretization  of 
the  Chandrasekhar  H-equation 


F(x)i  =  Xi  - 


N 


1  — -  y 

2  N  ^ 


3= 1 


HiXj 

Hi  +  Hj 


-l 


were  obtained  with  a  Newton-GMRES  solver  nsoli  from  [16].  We  set  the  relative 
and  absolute  tolerances  in  the  solver  to  10-12. 


c 

Fig.  3.1.  Bifurcation  diagram  of  the  H-equation  for  c  close  to  1.  The  parameter  c  is  equal  to 
0.9999179  at  the  point  marked  by  a  circle. 


Figure.  3.1  is  a  bifurcation  diagram  with  respect  to  the  parameter  c,  showing  a 
turning  point  at  c  =  1.  For  all  values  of  c  shown  in  this  continuation,  the  eigenvalues 
of  the  linearization  of  the  solution  are  close  to  1;  only  a  single  one  of  them  changes 
in  marked  way,  ranging  from  -0.9  on  the  upper  branch  through  zero  at  the  turning 
point  to  0.5  on  the  lower  branch. 

Figure.  3.2  shows  the  GMRES  performance  for  a  representative  value  of  c  (0.99991 
79),  close  to  the  turning  point,  marked  by  a  circle  in  Figure  3.1.  We  report  on 
computations  from  the  continuation  itself,  where  the  initial  iterate  was  the  standard 
linear  predictor,  and  from  a  second  stand-alone  computation  where  the  initial  iterate 
was  the  solution  plus  a  perturbation  function  0.05sin(:r).  In  Figure  3.2  we  plot  the 
convergence  history  for  Newton-GMRES  in  the  two  cases.  The  convergence  rates  and 
the  costs  of  the  solves  are  roughly  the  same. 

The  first  column  of  Table.  3.1  shows  the  ten  eigenvalues  farthest  from  1  for  the 
linearization  of  each  of  the  two  problems.  In  the  continuation  case,  there  is  one  more 
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Fig.  3.2.  Convergence  plot  for  the  steady  state  and  continuation  problem  of  H- equation 


eigenvalue  far  away  from  1  than  in  the  stand-alone  case.  In  both  cases  the  linearization 
clearly  fits  the  pattern  of  a  compact  perturbation  of  the  identity. 


Eigenvalues  of  steady  state  problem 

Eigenvalues  of  continuation  problem 

0.0207265 

-5.2022448 

0.9424114 

5.2008250 

0.9900391 

0.9746094 

0.9974043 

0.9828252 

0.9992541 

0.9977360 

0.9998164 

0.9991859 

0.9999612 

0.9998098 

0.9999926 

0.9999577 

0.9999987 

0.9999919 

0.9999998 

0.9999985 

Table  3.1 

Eigenvalues  for  the  linearized  steady  state  and  continuation  problems  of  the  H-equation  at 
c  =  0.9999179 


3.2.  The  Chafee-Infante  reaction  diffusion  problem.  We  then  solved  the 
steady  state  and  continuation  problems  for  a  discretization  of  a  dissipative  reaction- 
diffusion  PDE  in  one  dimension,  the  so-called  Chafee-Infante  problem, 

1  3 

‘U't  ttxx  u 


u  =  0,  x  £  [0, 7r] 
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with  boundary  conditions  u(0,t)  =  0,u(n,t)  =  0.  We  used  201  finite  difference 
discretization  points. 


Fig.  3.3.  Bifurcation  diagram  for  the  Chafee- Infante  reaction  diffusion  problem. 

Figure.  3.3  shows  the  bifurcation  diagram  for  this  problem  for  a  range  of  param- 
eter  values  (0  <  A  <  7)  where  up  to  five  different  spatially  structured  steady  states 
exist.  Our  computational  tests  were  performed  on  the  upper  stable  solution  branch 
close  to  A  =  2.  In  this  example  we  studied  both  the  steady  state  problem  (setting  the 
right-hand-side  of  the  finite  difference  equations  equal  to  zero)  and  the  time-stepper 
formulation  (our  integration  routine  was  odel5s).  The  steady  state  (equations)  /  fixed 
point  (time-stepper)  and  continuation  problems  were  solved  with  Newton-GMRES 
solver  nsoli  for  various  time  reporting  horizons.  The  results  for  A  =  2.1386697  are 
shown  in  Fig.  3.4.  Relative  and  absolute  tolerances  were  chosen  to  be  10~12.  The  ini¬ 
tial  guess  for  the  direct  solution  was  chosen  to  be  the  true  solution  plus  a  perturbation 
function  O.lsin(x)  . 

One  thing  should  be  made  clear  at  this  point.  Using  several  time  steps  of  an 
implicit  integrator,  with  the  concomitant  nonlinear  solves,  is  clearly  not  an  efficient 
way  of  solving  a  fixed  point  problem  (a  single  nonlinear  solve).  The  integrator  is  used 
here  as  a  “legacy  code” ,  a  code  that  one  cannot  modify.  It  is  in  the  context  of  such 
legacy  codes  that  our  approach  becomes  useful,  as  well  as  in  the  case  of  multiscale 
computations,  where  the  time  evolution  is  performed  by  a  simulator  at  a  different 
level  of  description  (e.g.  Lattice-Boltzmann  or  kinetic  Monte  Carlo). 

3.3.  Clustering  of  Eigenvalues.  Using  the  Finite  Difference  Method,  the  ODE 
system  from  discretizing  the  original  PDE  has  a  form  of  w(f)  =  f(u),u  £  Rn.  At  a 
steady  state  u*  ,  the  matrix  /„  has  n  eigenvalues:  <7j,  i  =  1, . . . ,  n.  If  we  solve  for  u* 
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Fig.  3.4.  Convergence  of  time- stepper  based  fixed  point  computation  for  the  steady  state  and 
the  continuation  of  the  discretized  RD  problem.  Every  circle  corresponds  to  a  Newton  step.  A  single 
“function  evaluation ”  consists  of  integration  over  one  time  reporting  horizon. 


with  a  time-stepper,  the  system  to  be  solved  is  w(0)  —  $t(w,  A)  =  0,  which  also  has  n 
eigenvalues:  1  —  eaiT ,  i  =  1, . . . ,  n. 

We  have  used  the  forward  in  time  time-stepper  to  compute  both  stable  and  un¬ 
stable  steady  states.  We  have  marked  the  unstable  steady  state  we  computed  using  a 
reporting  horizon  of  T  =  0.1  for  A  =  4.5710239  using  forward  in  time  time-stepping. 
When  the  parameter  A  is  equal  to  2.183867,  u*  is  a  stable  steady  state  and  all  the 
eigenvalues  are  negative.  When  we  increase  T,  all  the  eigenvalues  1  —  eaiT  are  ap¬ 
proaching  (“clustering  at”)  1.  Clustering  of  eigenvalues  is  known  to  be  beneficial  for 
GMRES  performance.  Conversely,  when  we  decrease  T,  eigenvalues  start  leaving  the 
cluster.  This  results  in  additional  GMRES  iterations. 

To  quantify  the  dependence  of  the  performance  of  the  iteration  on  T  we  use  the 
number  of  GMRES  iterations  at  the  final  step  in  Figure  3.5.  When  T  is  large,  we 
consistently  see  2  linear  iterations.  As  T  is  reduced,  we  see  an  increase  in  GMRES 
iterations,  as  expected.  In  this  particular  example,  when  T  =  1.78,  the  number  of 
GMRES  iterations  increases  to  3. 

In  at  attempt  to  quantify  this  further  we  identify  a  “cluster”  of  eigenvalues  that 
seems  to  be  correlated  with  the  performance  of  the  GMRES  iteration.  In  Fig.  3.6, 
we  treat  those  eigenvalues  in  the  interval  [0.84  1]  as  “in  the  cluster”.  The  number 
of  eigenvalues  outside  the  cluster  (or  smaller  than  0.84)  and  the  number  of  function 
evaluations  needed  to  finish  the  last  Newton  step  are  compared  in  Tab.  3.2  below.  A 
clear,  strong  correlation  emerges. 

3.4.  Conclusion.  We  have  extended  and  sharpened  results  for  the  performance 
of  GMRES  for  discretizations  of  compact  fixed  point  problems  to  the  special  class 
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Fig.  3.5.  Convergence  of  time- stepper  based  fixed  point  Newton-GMRES  computation  for  dif¬ 
ferent  time  reporting  horizons.  Every  circle  corresponds  to  a  Newton  step,  and  m  is  the  number  of 
function  evaluations  for  the  last  Newton  step. 


Fig.  3.6.  Twenty  smallest  eigenvalues  of  the  linearized  time- stepper  at  the  steady  state  for 
different  time  reporting  horizons;  the  dashed  line  corresponds  to  the  eigenvalue  that  first  leaves  the 
cluster  when  T  =  1.78  . 
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T 

4 

2 

1 

0.5 

0.3 

0.1 

0.07 

0.04 

0.02 

Number  of 
eigenvalues  outside 
the  cluster 

0 

0 

1(2)* 

2 

3 

6 

7 

10 

14 

Function  evaluations 
needed  for  clustering 
eigenvalues 

2 

2 

2 

2 

2 

2 

2 

2 

2 

Actual 

total  function 
evaluations 

2 

2 

5 

6 

7 

11 

8 

12 

16 

Table  3.2 


of  problems  that  arise  in  a  time-stepper  context.  The  key  feature  leading  to  the 
enhancement  of  GMRES  performance  is  the  compactification  of  the  spectrum.  In 
previous  work,  a  Green’s  function  based  reformulation  of  the  steady  state  problem 
for  elliptic  PDEs  [7]  led  to  a  linearization  that  was  a  compact  perturbation  of  the 
identity,  and  an  efficient  solution  via  Newton-GMRES.  Time-steppers  compactify  the 
spectrum  of  the  linearization  in  a  natural  way,  and  their  properties  can  be  exploited  to 
obtain  accurate  bounds  on  the  convergence  rates  of  the  linear  iterations  in  a  Newton- 
GMRES  continuation.  We  show  that  the  additional  equation  in  a  pseudo-arclength 
formulation  of  a  parameter-dependent  family  of  nonlinear  equations  adds  at  most  two 
GMRES  iterations  when  the  eigenvalues  are  well  separated,  and  obtain  a  bound  on 
the  convergence  in  terms  of  the  separation  of  the  spectrum  and  the  dimension  of  the 
slow  subspace  (associated  with  a  slow/inertial  manifold  for  the  dynamics). 

We  reported  on  numerical  results  that  support  the  theory,  first  with  an  integral 
equation,  for  which  we  can  numerically  demonstrate  that  all  but  two  of  the  eigen¬ 
values  of  the  linearization  lie  in  a  tight  cluster  about  1.  The  second  example  was 
a  parametric  study  of  the  steady  state  of  a  discretized  parabolic  partial  differential 
equation  implemented  through  time-steppers. 

The  “natural”  compactification  of  the  spectrum  using  time-steppers  provides  a 
natural  connection  with  the  performance  of  matrix-free  iterative  methods.  This  com¬ 
pactification  may  prove  useful  in  writing  computational  wrappers,  that  will  accelerate 
the  convergence  of  legacy  dynamic  simulators  to  stationary  states.  We  also  expect  this 
compactification  to  assist  in  writing  computational  wrappers  that  will  assist  dynamic 
simulators  at  a  different  level  of  model  description  (e.g.  kinetic  Monte  Carlo,  Brow¬ 
nian  Dynamics  or  Molecular  Dynamics  codes  [11,20,24])  to  locate  stationary  states 
and  perform  continuation/bifurcation  analysis  of  macroscopic  system  observables  in 
the  so-called  equation-free  framework  [19]. 
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